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SUMMARY 

Full-field  strain  measurement  techniques  are  based  on  computing  the  spatial  derivatives  of  numerical  or 
functional  approximations  of  the  underlying  displacement  fields  extracted  from  digital  imaging  methods. 
These  methods  implicitly  assume  that  the  medium  satisfies  the  strain  compatibility  conditions,  which  are 
only  true  in  the  case  of  a  continuum  body  that  remains  continuum  throughout  its  deformation  history.  In  the 
present  work,  we  introduce  a  method  that  can  be  used  to  calculate  the  strain  components  directly  from  typ¬ 
ical  digital  imaging  data,  without  the  need  of  the  continuum  hypothesis  and  the  need  for  displacement  field 
differentiation.  Thus,  it  enables  the  measurement  of  strain  fields  from  imaged  surfaces  that  may  or  may  not 
contain  discontinuities.  Numerical  comparisons  are  performed  on  the  basis  synthetic  data  produced  from  an 
analytical  solution  for  an  elastically  orthotropic  open-hole  domain  in  tension.  For  performance  comparison 
purposes,  the  mean  absolute  error  distributions  are  calculated  for  the  cases  of  both  the  traditional  meshless 
random  grid  method,  and  the  direct  strain  method  introduced  herein.  It  is  established  that  the  more  refined 
representation  of  strain  provided  by  our  present  approach  is  more  accurate  everywhere  in  the  domain,  but 
most  importantly,  near  its  boundaries.  Published  2013.  This  article  is  a  US  Government  work  and  is  in  the 
public  domain  in  the  USA. 
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1.  INTRODUCTION 

Calculation  of  strain  by  differentiation  of  displacement  fields  is  a  typical  practice  in  full-field  mea¬ 
surement  methods  [1-15].  In  addition,  there  are  techniques,  such  as  Shearography  [16-18]  and 
Moire  interferometry  [19,20]  that  exploit  implicit  differentiations  of  the  displacement  fields,  com¬ 
pute  only  a  subset  of  the  strain  components,  and  require  specialized  instrumentation,  are  impractical 
for  nonplanar  cases  and  are  sensitive  to  out  of  plane  motions. 

Strain  field  estimation  based  on  differentiation  of  displacement  fields  utilized  by  these  techniques 
implies  that  the  strain  compatibility  equations  remain  satisfied  for  the  entire  domain  of  observation 
over  the  entire  loading  history  such  as  the  continuum  hypothesis  remains  valid.  However,  when 
surface  discontinuities  or  strain  localization  occur  because  of  damage  initiation,  the  continuum 
hypothesis  is  no  longer  valid  throughout  the  entire  field  of  observation.  Therefore,  full-field  mea¬ 
surement  results  based  on  the  anticipated  validity  of  the  continuous  hypothesis  yield  both  inaccurate 
and  indeterminate  measurements,  thus  leading  to  false  qualitative  and  quantitative  conclusions. 

Our  experience  with  the  utilization  of  high-throughput  multidegree  of  freedom  automated  mecha- 
tronic  testing  machines,  which  are  capable  of  loading  specimens  multiaxially  in  conjunction  with 
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energy-based  inverse  constitutive  characterization  methodologies  for  composite  materials  [21-26], 
has  underlined  the  need  for  both  automated  full-field  data  processing,  but  also  of  an  accu¬ 
rate  encapsulation  of  deformation  fields.  Such  concerns  become  increasingly  more  important 
in  regions  of  the  specimen  that  are  close  to  or  have  undergone  failure,  and  accordingly,  the 
medium  can  no  longer  be  considered  as  a  continuum.  In  these  areas,  the  microcracks  and  cracks 
developed  break  the  validity  of  the  continuum  hypothesis  and  the  associated  strain  compatibility 
relations. 

Another  important  issue  arising  from  the  derivation  of  strain  from  a  full-field  displacement  repre¬ 
sentation  is  that  of  relative  lower  accuracy  near  and  on  the  boundaries  of  the  implemented  analytical 
representation  [13,27],  As  the  strain  components  are  defined  in  terms  of  the  derivatives  of  the  dis¬ 
placement  components,  it  is  expected  that  the  analytical  approximation  of  the  strain  components 
will  be  somewhat  problematic,  as  the  displacement  approximation  does  not  impose  spatially  driven 
constraints  on  the  formulation  of  its  derivatives.  Although  it  is  possible  to  reformulate  the  meshfree 
approximation  of  the  displacement  fields  to  take  into  consideration  such  information,  this  would  be 
impractical  and  would  require  considerable  and  customized  effort  to  manually  guide  the  approxi¬ 
mation  from  an  algorithmic  perspective  in  a  manner  that  is  aware  of  the  boundaries,  and  it  remains 
general  without  need  for  customization  from  geometry  to  geometry. 

The  final  issue  of  concern  arises  from  the  fact  that  if  there  is  noise  in  the  displacement  fields,  their 
differentiation  will  only  amplify  it  for  the  evaluation  of  the  strain  fields. 

The  challenging  character  of  these  issues  has  driven  us  to  consider  a  new  approach  in  calculat¬ 
ing  the  full  field  of  strain  quantities  from  full-field  digital  images  of  deforming  specimens.  The 
numerical  and  analytical  method  presented  in  this  paper  succeeds  in  representing  the  strain  ten¬ 
sor  field  directly  from  directional  engineering  strain  quantities  without  requiring  or  enforcing  the 
satisfaction  of  the  compatibility  conditions. 

Because  this  paper  aims  at  introducing  direct  strain  imaging  (DSI)  as  a  new  full-field  mea¬ 
surement  method,  we  first  describe  the  steps  involved  with  the  implementation  of  the  traditional 
meshless  random  grid  (MRG)  method  [5-15].  This  will  help  delineate  the  differences  and  provide 
the  basis  of  performance  comparison  between  MRG  and  DSI  methods. 

The  typical  experimental  procedure  for  measuring  the  full  field  of  deformation  quantities  accord¬ 
ing  to  most  full  field  methods  in  general,  and  MRG  method  in  particular,  can  be  outlined 
as  follows: 

1 .  A  specimen  is  marked  with  an  appropriate  visible  pattern  that  consists  of  a  random  distribution 
of  dots  distinguishable  from  the  background. 

2.  If  the  experiment  is  to  take  into  consideration  out-of-plane  motion,  two  or  more  cameras  are 
used  so  that  the  deformation  is  stereoscopically  reconstructed.  The  projective  characteristics 
of  the  cameras  are  identified  through  an  appropriate  calibration  procedure  [8,28] 

3.  The  specimen  is  placed  in  the  mechanical  testing  machine,  and  one  image  per  camera  is 
captured  in  the  undeformed  configuration  prior  to  the  initiation  of  the  loading  sequence. 

4.  While  the  experiment  is  taking  place,  successive  images  of  the  deforming  specimen  are 
captured. 

5.  The  images  are  processed  and  for  each  frame,  the  coordinates  of  appropriate  points  (nodes) 
are  calculated.  Those  nodes  may  be,  for  example,  geometric  centroids  of  dots  (for  the  case 
of  grid  methods),  appropriate  boundaries  of  geometric  entities,  correlated  subregions  (for  the 
case  of  Digital  Image  Correlation  (DIC)),  and  so  on. 

6.  Using  an  interpolation  or  approximation  scheme,  a  representation  of  the  displacement 
field  is  obtained.  In  the  MRG  method,  the  field  is  represented  by  a  continuous  meshless 
approximation. 

7.  The  strain  at  any  point  in  the  domain  is  obtained  by  differentiation  of  the  displacement  fields 
based  on  the  definition  of  the  strain  measure  as  a  function  of  the  displacements. 

Noise  introduced  by  various  sources  is  the  most  dominant  source  of  error  and  plagues  all  full- 
field  measurement  methods.  In  its  presence,  the  MRG  method  has  been  shown  to  perform  very 
well  [13]  compared  with  other  methods.  The  main  reason  for  its  improved  performance  is  that  the 
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approximation  scheme  of  the  meshless  representation  works  naturally  as  a  filter  that  deals  excep¬ 
tionally  well  with  the  image  acquisition  noise.  Unfortunately,  because  not  all  noise  can  be  removed, 
there  is  still  room  for  improving  the  accuracy  levels  of  full-field  techniques  for  certain  experimental 
configurations. 

Following  this  introduction,  the  paper  continuous  with  a  brief  description  of  the  overall  problem 
at  hand,  and  introduces  the  DSI  method  by  describing  a  meshless  approximation  algorithm  applied 
on  strain  tensor  quantities  derived  from  distinct  engineering  strain  component  measurements.  A 
proof  of  the  consistency  of  this  approximation  follows.  The  paper  continues  with  presentation  of 
synthetic  experiments  that  are  developed  to  explore  the  performance  of  the  method  on  well  estab¬ 
lished  full-field  representations,  based  on  the  deformation  of  an  ortho  tropic  plate  with  and  open  hole 
in  tension.  The  paper  closes  with  conclusions  of  the  current  work  and  future  plans. 


2.  PROBLEM  STATEMENT 

The  engineering  strain  between  pairs  of  adjacent  centroids  of  dots  (nodes)  along  the  directions 
defined  by  them,  will  be  the  digitally  determined  data  input  for  the  proposed  analysis.  A  domain  G 
in  the  undeformed  configuration  is  shown  in  Figure  1(a),  populated  with  a  number  of  such  nodes. 
The  same  domain  in  the  deformed  configuration  is  shown  in  Figure  1(b).  For  a  set  of  n  nodes  in  the 
vicinity  of  an  interest  point  described  by  a  position  vector  w  =  {xw,yw}T ,  it  is  possible  to  form 
a  number  of  engineering  strain  quantities  e,-,-  =  (l^  —  /,-,)/ lij,i  —  1  ...  n  —  1,  j  —2  <  j 

[29]  in  terms  of  the  distances  /y ,  /(- ■  in  the  initial  and  deformed  configurations,  respectively.  Each 
of  those  quantities  can  be  considered  as  the  average  of  the  strain  component  in  the  direction  of  the 
line  connecting  the  two  nodes  over  the  respective  line  segment.  Given  those  engineering  strains  and 
the  coordinates  of  the  nodes,  we  seek  to  calculate  the  strain  tensor  at  any  point  w  in  the  domain.  It 
is  important  to  note  that  similar  to  traditional  meshless  methods  [30-32]  the  vicinity  of  a  point  of 
interest  w  =  {xw ,  yw  }T  is  bound  by  an  appropriate  domain  Of  support  (DOS),  and  in  this  work,  we 
will  consider  the  DOS  to  be  circular. 

The  solution  of  this  problem  might  be  at  first  considered  as  the  extension  of  the  typical  three  gauge 
rosette  strain  calculation  [29].  Flowever,  this  would  only  hold  true  if  the  strain  field  is  considered 
constant  throughout  the  region  of  interest.  Indeed,  in  calculating  the  strain  from  a  strain  rosette,  it 
is  assumed  that  the  engineering  strains  measured  by  each  directional  strain  gauge  of  the  rosette,  all 
occur  at  the  same  point.  Although  this  might  be  true  for  certain  strain  gauge  rosette  configurations, 
it  is  certainly  not  the  case  under  the  requirements  of  the  problem  we  are  attempting  to  solve  because 
as  is  evident  from  Figure  1 ,  the  engineering  strain  quantities  cannot  be  considered  as  been  measured 
at  the  same  point. 


(a)  Undeformed  (b)  Deformed 

Figure  1 .  Deformation  configurations  and  distances  of  selected  node  pairs. 
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It  should  be  noted  that  although  the  method  presented  herein  is  used  in  the  context  of  full-field 
strain  measurements  using  digital  imaging  information  in  the  visual  spectrum,  it  can  be  also  applied 
to  any  context  where  appropriate  data  can  be  acquired.  For  example,  engineering  strain  data  may 
be  obtained  by  quantifying  separation  distances  between  neighboring  markers  via  the  use  of  ultra¬ 
sound  or  X-ray  tomography  techniques,  and  the  respective  full-field  strain  tensor  could  be  calculated 
between  appropriate  volumetric  speckles  by  the  method  presented  in  this  paper.  The  same  is  true  for 
all  methods  that  can  identify  marker  separation  on  surfaces  of  deformable  bodies. 

In  the  following  section,  we  will  present  a  solution  to  this  problem  by  the  means  of  a  spatial 
moving  least  squares  approximation. 

3.  DIRECT  STRAIN  APPROXIMATION 


Our  aim  is  to  identify  the  strain  tensor  components  at  any  point  described  by  a  position  vec¬ 
tor  w  =  {xw,  yw  }  T  within  a  deformable  domain.  An  approximation  of  the  components  of  the 
two-dimensional  strain  tensor  at  any  point  x  =  {x,  y}T  in  the  vicinity  of  w  can  be  expressed  by 


m 

ejx(w;  x)  =  ^2  Pj(x)a*x( w)  =  pT(x)ax*(w) 
j 

m 

4(w;x)  =  £>,  (x)af  (w)  =  pT(x)a^(w)  .  , 
j 

m 

s^y(w:x)  =  ^Pj{x)aY  (w)  =  pT(x)axy(vf ) 
j 


(1) 


where  p(x)  is  a  vector  of  basis  functions  that  can  be  chosen  to  consist  of  in  monomials  of  the  lowest 
orders  to  ensure  minimum  completeness.  A  polynomial  basis  of  order  M  <  m  has  the  general  form 

P(x)  =  p(x,y)  =  {l,x,y,xy,x2,y2,...,xM,yM}T  .  (2) 

Such  a  polynomial  basis  can  be  constructed  by  concatenation  of  the  complete  order  terms  using  the 
Pascal  triangle  of  monomials  [30, 32].  It  should  be  noted  that  the  polynomial  basis  could  be  chosen 
to  be  different  for  each  of  the  strain  components,  but  this  would  add  in  algebraic  complexity  that  is 
not  presently  justified.  In  addition,  a“  (w)  and  a  =  xx,  yy,xy  are  vectors  of  coefficients  specific  to 
every  point  w  of  the  domain, 

a“  (w)  =  {a% ,  a “ , . . . ,  a“  }T  ,  a  =  xx,  yy ,  xy .  (3) 


In  contrast  to  the  MRG  method  [8],  here  we  are  seeking  to  directly  calculate  the  strain  com¬ 
ponents.  Therefore,  we  first  need  to  identify  the  engineering  strains  between  a  number  of  node 
combinations. 

The  normal  strain  at  an  angle  9  relative  to  the  global  coordinate  system  is  given  by  the  well-known 
transformation, 

sxx  +  £yy  £xx  ~  £yy  ,,, 

exx'  = -  H - cos  26  +  Sxy  Sin  29 .  (4) 

Using  Equations  (4)  and  (1),  we  can  derive  the  function  that  approximates  the  xx  strain  tensor 
component  at  a  point  x  =  {x,  y}T  on  a  coordinate  system  Ax' y'  forming  an  angle  9  with  the  global 
coordinate  system  Axy  (Figure  2)  in  the  vicinity  of  w,  as  follows: 

h  ,  ,  pT(x)ax*(w)  +  pT(x)a-v>'(w) 

e"x,(w:x)  = - - - + 

pT(x)axx(w)-pT(x)a^(w)  (5) 

H - cos  29+ 

2 

+  pT(x)ax;v  (w)  sin  29, 


Published  2013.  This  article  is  a  US  Government  work  and 
is  in  the  public  domain  in  the  USA. 


Int.  J.  Numer.  Meth.  Engng  2013;  95:313-330 
DOI:  10.1002/nme 


DIRECT  STRAIN  TENSOR  APPROXIMATION 


317 


Figure  2.  Strain  component  sx  ,  of  the  approximation  in  the  direction  of  AB  relative  to  the  local  coordinate 

system  at  point  A. 


or 


sxx,( w;  x)  =  pT(x)axx(w)  cos2  9+ 

+  pT  (x)ayy  (w)  sin2  9  +  pT(x)axy(w)  sin  29. 

Within  the  domain  Q  (Figure  2)  and  in  the  vicinity  of  a  point  expressed  by  the  vector  w,  we 
assume  the  existence  of  two  neighboring  nodes  A  and  B  expressed  by  the  vectors  x^  =  {xA,  }’a  }t 
and  xB  —  {rs,  yB  }t,  respectively.  The  mean  strain  between  these  two  nodes  as  derived  by 
Equation  (6)  can  be  defined  by 


£x*'(W:X^’XB)  = 


fo  e*x,(w;s(xA,xB;X )) 


9s(x^,x^;A) 


3A 


dA 


|xs  -x^ 


(7) 


where  s(x^,  x^ ;  A)  is  the  parametric  expression  of  the  line  segment  passing  from  both  A  and  B  and 
is  given  by 


s(x^,xB; A)  =  (1  -  A)x^  +  XxB  — 

—  {(1  —  +  XxB,  (1  —  A)y^  +  Ays},  A  e  [0, 1], 


(8) 


It  is  easy  to  prove  that 


ds(xA,xB:X) 

3A 


|xg  -xA\, 


(9) 


and  hence  Equation  (7)  simplifies  to 

shxx,(w,xA,xB)  =  f  sxx^w.s(xA,xB-.X))dX.  (10) 

Jo 

The  line  integral  in  Equation  (10)  can  then  be  derived  by  combining  Equations  (6)  and  (8), 

/  sxx,(yr.  s  (xA,xB\X))  dA  = 

Jo 


—  [  [pT(s)axx(w)  cos2  Od- 
Jo 

+  pT(s)a;>;>’(w)  sin2  9  +  p^sja^w)  sin  26]  dA. 


(ID 


Published  2013.  This  article  is  a  US  Government  work  and 
is  in  the  public  domain  in  the  USA. 


Int.  J.  Numer.  Meth.  Engng  2013;  95:313-330 
DOI:  10.1002/nme 


318 


A.  ILIOPOULOS  AND  J.  G.  MICHOPOULOS 


The  angle  6  can  be  calculated  by  the  nonambiguous  arctangent  of  the  vector  xB  —  xA .  Splitting  up 
the  integral  and  distributing  the  independent  terms  outside  the  integrands  yields 


/' 


e^x,(w;s(x^,xg;A))dA  = 


cos2  9  /  pT(s  (x^x#;  A))  dA 


’/:» 

•L 


sin2  6  /  pT(s  (x^xg;  A))  dA 


sin  26 


/v 


(s  (xA,xB:  A))  dA 


axx(w)+ 
a^(w)+ 
axy  (w). 


(12) 


The  vector  integrals  in  Equation  (12)  can  be  easily  calculated  for  a  given  polynomial  basis  and 
can  be  derived  from  the  terms  of  Table  I.  If  a  full-order  polynomial  is  chosen,  the  integral  basis 
vector  can  be  constructed  by  concatenation  of  the  respective  terms.  For  example,  if  a  six-term  basis 
vector  is  chosen  to  be  p(x)  =  p(x,  y)  =  { 1,  x,  y,  xy,  x2,  y2}T,  the  integral  basis  vector  will  be 


=  f 


r  T(xA,xB)  =  pT  (s(x^,xB;  A))dA 


=  1 1,  ]-(xA  +  xB),  ~ {yA  +  Vb),  ^  (xA  (2 yA  +  yB)  +  xB(yA  +  2 yB)) ,  (13) 


+  xAxB  +  xB),  ^(yA  +  y Ay b  +  y\)  \  ■ 


Table  I.  Integral  basis  monomials. 


Complete 

order 

Number  of 
terms 

Additional  terms 

Constant 

1 

1 

Linear 

3 

\k(xA  +xB), \{yA  +  Tb)} 

Bilinear 

4 

\(XA  (2yA  +  yB)  +  xB(yA  +  2.vg» 

Quadratic 

6 

{ 3  (xaxb  +  XA  +  X%),\  {yAvB  +  y\  +  y% ) } 

Quadratic 

8 

{ T2 (XA  (3TA  +  yB )  +  2xAxb (yA  +  yB)  +  X2  (yA  +  3 yB )) , 

T2  (yA  (3xa  +  xB)  +  2 yAyB(xA  +  xB)  +  y\(xA  +  3xs))| 

Cubic 

10 

{\(xA  +  xB)  (x\  +  X2B) ,  4 (yA  +  yB)  (yA  +  y|)} 

Quartic 

15 

{  Jo  (XA  (4-v4  +  }'B )  +  x2axb  (3 yA  +  2 VB )  +  xAx2  (2 yA  +  3 yB )  + 
+xB(yA  +  4yB)) , 

30  (xA  {3yA\’B  +  6 yA  +  y%)  +  xAxB  (4 yAyB  +  3 y\  +  3y|)  + 
+X2B  (iyAyB  +  yA  +  6y|)) , 

Jo  (yA  (4xa  +  xB)  +  y\yB( 3xa  +  2xB)  +  yAyB(2xA  +  3xB)  + 
+yB(xA  +  4xB)), 

j  (xAXB  +  x2ax2b  +  XAXB  +  XA  +  XB) , 

5  (y\ys  +  y\y\  +  vavI  +  yA  +  yB)} 
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The  choice  of  an  appropriate  polynomial  basis  will  yield  an  integral  vector  basis  rT(x^,  xg).  Thus, 
Equation  (12)  can  be  written  as 


or 


e*x'(w:  s  (xA,xg;  A))  = 

=  [r'(x/1,xB)cos2  0]  ax*(w)+ 

+  [rT(x^,xg)  sin2  0]  aJ’;>;(w)+ 
+  [rT(x^,xB)sin20]aX3;(w), 


(14) 


eJx,(w;s(x^,xB;A))dA  = 

=  qL(xA,xg)axx(w)  +  qJ>,(x4,xg)aJ,>’(w)  + 
+  q  Jj,  (x^ ,  xg  )a*-)'  (w) , 


(15) 


with 


qL(X/f,xg)  =  rT(xA,xg)cos2  01 

qyj(xA,xg)  =  rT(x^,xg)  sin2  6  l.  (16) 

qJy(xT,xg)  =  rT(x^,xg)  sin 20  J 


Collecting  the  terms  of  Equation  (15)  into  larger  vectors  and  substituting  in  (7)  yields 

8*  ,(w;x^,xg)  =  qT(x^,xg)a(w), 


(17) 


with 


qT  (XA ,  xg )  =  {qkB}  = 

=  {qJx(xA,xg),  q^(x^,xg),  qjy(x^,xg)} 

/  ,  T  >  .  (IB) 

a(w)  =  {ak}  =  |(axx  (w))T,  (a^  (w))T  ,  (a**  (w))TJ 

k  —  1 . . .  3m 

Equation  (17)  describes  the  approximated  mean  strain  between  any  two  nodes  in  the  domain  at 
the  vicinity  of  point  w.  When  the  coordinates  of  the  nodes  are  determined  using  a  digital  imaging 
technique,  it  is  possible  to  calculate  the  real  mean  strain  (aka  the  engineering  strain)  between  any 
two  of  the  n  nodes  as 


C/  —  t'xx'  (X;  ’  x  / ) 


I'ii  -  hj 


li 


^-1  = 


— - J—  —  l,z  —  1 ...  n  —  1,  j  =  2 ...  n,  i  <  j , 

x7  xy  I 


(19) 


where  x,-  and  x;  are  the  vectors  representing  the  positions  of  nodes  i,  j  at  the  undeformed  configu¬ 
ration,  and  x(  and  x  ■  are  the  same  vectors  in  the  deformed  configuration.  The  scalar  quantities  etJ 
should  not  be  thought  of  as  the  elements  of  a  matrix,  but  rather  as  a  set  of  n(n  —  l)/2  numbers. 
Although,  not  of  interest  in  the  current  analysis,  the  latter  can  be  thought  as  functions  of  x,  and  x7 
and  time.  The  condition  expressing  that  i  should  be  greater  than  j  ensures  the  requirement  that 
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no  combination  is  considered  twice.  A  functional  of  weighted  residuals  can  now  be  constructed  by 
using  Equations  (17)  and  (19),  as  follows: 


n—l,n 

J  —  E  W^w  — 

i=l,j=2,i<j 
n—\,n 

1  w  - 


X;  +  X; 


)(4'  (w;xi,x7)-e^7(xi,xy))' 


(20) 


E  w{' 

--l,j=2,i<j  V 


X;  +  Xj  >  ,  T 


(*r  (x; , x7)  a(w)  —  eij)  , 


where  W(yv—  (x,-  +xy)/2)  =  Wij  ^  0  is  a  weight  function  that  decreases  with  distance  as  introduced 
by  [30, 32].  A  common  weight  function  is  the  cubic  spline  weight  function  given  by, 


W(  d)  =  W(d)  = 


- Ad2  +  4d3 

3 

4  a  A  .  4  A-  |  /V 

-~4d  +  4d2--d3  ,-<dH  1 
3  3  2 


(21) 


0 


,  d  >  1 


with  d  —  -j-,  the  normalized  distance,  and  dw  the  smoothing  length  that  is  usually  equal  to  the 
extent  of  the  domain  of  support. 

It  should  be  noted  that  any  strain  component  could  be  used  in  the  functional  of  Equation  (20)  and 
that  all  of  them  are  present  in  this  equation  because  of  Equation  (4).  The  notation  i  —  1 J  =2,2  <  j 
represents  the  set  of  all  combinations  between  the  indexes  i  and  j .  The  least  squares  minimization 
of  the  residual  defined  by  Equation  (20)  will  provide  the  values  of  coefficients  in  the  vector  a(w)  so 
that  shxx,  (w;  x, ,  xy )  optimally  approximates  sxx>  (x(- ,  x7  )  .  To  minimize  the  residual  in  Equation  (20) 
with  respect  to  a,  we  form  the  derivatives 


9  J  dJ 

—  =0  =*>  - —  =  0,  k  =  1 . . .  3m 
9a  da/c 


(22) 


that  can  be  expanded  to 

3J  3  Etlj, ix,  Wi  (£"■  <#«'  -  ‘u) 
dak 


dak 

—  2  Wij  I  Y'  qlrJar  I 


n— l,n 


\r=  1 
/  3  m 


dak 


i  =  l,j=2,i<j 
n — l,n 

=  2  W‘J  ( J2  Ur  -  ea )  vlk  = 

i  =  l,j=2,i<j  \r= 1  / 

n—l,n  3  m  n—l,n 

=  2  E  WaJ2(lk<lriar-2  E  W‘jdkeV- 

i=l,j=2,i<j  r= 1  i=l,j=2,i<j 

Introducing  this  expression  to  Equation  (22)  yields 

n—l,n  /  3  m  \  n—l,n 

E  (  Wl>  E  dr’arj  =  E  W‘i 9lk  eH  ■ 

i  =  l,j  =2,i < j  \  r=  1  /  i  =  \J =2,i < j 

In  matrix  notation,  the  last  expression  can  be  written  as 

Q(w)a(w)  =  B(w)ew, 
or 


(23) 


(24) 


(25) 


a(w)  =  Q  1(w)B(w)ew  =  M(w)ew, 


(26) 
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where 


ew 


{^21,  P31,  ^32,  P41.  ^42,  •  •  •  ,  Pn(n-l)} 


T 


(27) 


is  a  vector  that  collects  all  the  measured  engineering  strains  in  the  vicinity  of  w,  B  is  a  matrix  that 
has  the  form 


B(w)  —  [b21,b3i,b32,b4i,b42,  .  .  .  ,b„(„_i)] 
by(w)  =  W  (w-  q(x,- ,xj) 

and  Q(w)  is  defined  by 

n-l.n  ,  x 

Q(W)=  wU-  ^  )  q(x« , x7 )qT (x, ,  x7 ) .  (29) 

i  =  l,j=2,i<j  x  J 


For  the  system  of  Equation  (25)  to  be  fully  determined  the  number  of  equations  has  to  be 
equal  to  3 m,  where  m  is  the  number  of  monomials  in  the  polynomial  basis.  For  most  cases, 
the  strain  values  are  to  be  acquired  by  all  the  combinations  between  n  nodes,  and  therefore: 

yy(yi 1)  772  ^  1  1  <■  /  \ 

v  2  ’  —  3 m  - >  n  —  2  ( 1  +  v  1  +  24m  ) .  However,  an  overdetermined  system  would  lead 

to  noise  reduction,  and  thus  is  preferable.  This  can  be  expressed  by 

-(l  +  Vl  +  24 m'j  .  (30) 

This  inequality  can  be  used  in  conjunction  with  additional  considerations  for  the  identification  of 
the  extent  of  the  domain  of  support  and  is  a  topic  of  a  follow-up  study. 

The  matrix  M(w)  =  Q-1(w)B(w)  will  be  of  size  3m  x  n,  where  each  row  corresponds  to  an 
element  of  the  vectors  described  in  Equation  (3). 

Partitioning  this  matrix  vertically  in  three  submatrices  of  dimensions  m  xn  yields 


M(w)  = 


Mxx(w) 
M^(w) 
Mx  j  (w) 


It  is  now  possible  to  express  the  coefficients  vectors  as 


axx(w)  =  Mxx(w)ew  j 
ayy(w)  =  \Iv>,(w)ew  1 
axy  (w)  =  MXJ,  (w)ew  J 


(31) 


(32) 


Substitution  of  these  expressions  in  Equation  (1)  yields  the  final  form  of  the  fully  defined  strain  field 
approximation  expressions 

eL(w)  =  ehxx(y/\y/)  =  <t>xx(w)ewj 

£hyy(y/)  =  Syy  (w:  w)  =  Qyy  (w)Cw  l  ,  (33) 

sxy(w)  =  sxy  (w:  w)  =  4»x>,(w)ew  J 

where  $«(w)  =  pT(w)Ma(w),  are  shape  functions  of  the  strain  components  for  the  point  w  for 
a  =  xx,  yy,  xy. 


4.  CONSISTENCY  TEST 

A  desired  but  not  necessary  property  of  approximations  is  that  of  consistency.  That  is  to  say,  given  a 
field  that  is  expressed  by  a  specific  basis  function,  then  it  should  be  reproducible  in  an  exact  manner 
by  the  approximation  when  using  this  basis  function. 
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Following  a  similar  approach  as  in  [30],  we  consider  that  the  actual  fields  of  the  components  of 
the  strain  tensor  of  a  domain  are  given  by 


-yy 


°xy 


(*)  = 

j 

k 

(X)  =  J2  Pj(x)cy/(x) 


(x)  =  ^Pj(x)cXjy(x) 


(34) 


We  also  assume  that  the  same  fields  are  to  be  approximated  by  functions  of  the  form  given  by 
Equation  (1)  with  in  —  k  and  the  same  basis  with  Equation  (34),  that  is, 


xx(x)  =  J2PjWaTW 
j 
k 

h  (x)  =  £>,  (x)af(x) 
j 
k 

(X)  =  Y^Pj{x)af  {x) 


Jyy 


-xy 


(35) 


It  is  obvious  that  J  in  Equation  (20)  will  be  equal  to  0  for  a“(w)  =  c“(x)  and  because  J  is  pos¬ 
itive  definite,  this  value  will  represent  the  minimum  minimorum.  In  addition,  because  the  sum  of 
Equation  (35)  yields  a  linear  combination  of  the  coefficients;  this  minimum  will  also  be  unique. 


Hence,  replacing  a“  with  c“  in  (35)  yields 


eL(x)  =  £xx(x),  £hxx(x)  =  £yy  (x),  Ehxx(x)  =  £Xy(x). 


(36) 


This  proves  that  the  approximated  fields  of  Equation  (35)  and  hence  Equation  (1)  reproduce  exactly 
the  fields  described  by  Equation  (34). 


5.  CONDITION  NUMBER 


As  indicated  in  [33],  the  effects  of  the  condition  number  in  meshless  surface  fitting  are  critical  only 
in  the  sense  of  invertibility  and  do  not  have  a  proven  relationship  with  the  effectiveness  of  the  regres¬ 
sion.  Accordingly,  in  the  case  of  the  DSI  method,  the  uniqueness  of  the  solution  -  as  expressed  by 
the  condition  number  -  is  not  of  utmost  importance  either.  However,  it  is  worthwhile  mentioning 
that  the  condition  number  should  not  be  outrageously  high  to  allow  the  approximation  of  the  tensor 
surface  from  a  numerical  perspective.  To  demonstrate  this  point,  we  consider  the  system  matrix  of 
interest  as  indicated  in  Equation  (25)  is  Q. 

To  study  the  effects  of  the  domain  of  support  and  the  number  of  monomial  terms  used  in  the 
polynomial  basis,  a  number  of  numerical  tests  were  performed.  A  random  distribution  of  nodes  was 
generated  in  a  square  domain  of  a  unit  length  as  shown  in  Figure  3.  The  matrix  Q  was  calculated 
using  Equation  (29)  for  a  number  of  domain  of  support  radius  values  and  for  a  number  of  monomial 
terms  of  the  polynomial  basis  function  of  Equation  (2). 

The  condition  number  of  matrix  Q  is  plotted  versus  the  domain  of  support  radius  in  Figure  4.  In 
this  plot,  it  is  observed  that  for  polynomial  basis  with  one  and  three  terms,  the  condition  number 
is  very  low,  indicating  very  stable  solutions.  For  higher  number  of  monomial  terms,  the  respective 
condition  numbers  tends  to  be  larger  but  in  acceptable  ranges  especially  when  compared  with  other 
meshless  methods  [33]. 
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Figure  3.  Randomly  generated  points  in  a  1  x  1  domain  used  in  the  study  of  the  condition  number  of  Q  of 

Equation  (25). 


Figure  4.  Condition  number  of  matrix  Q  versus  radius  of  domain  of  support  for  polynomials  with  various 

numbers  of  monomial  terms. 


6.  SYNTHETIC  EXPERIMENTS 

Because  DSI  is  a  consistent  approximation,  it  is  expected  that  it  can  exactly  reproduce  simple  poly¬ 
nomial  fields  and  hence,  it  is  more  interesting  to  validate  it  in  more  realistic  and  complex  situations. 
To  this  end,  a  set  of  synthetic  validation  experiments  was  conducted  on  a  well-known  deformation 
field  resulting  front  the  tensile  loading  of  an  orthotropic  body  with  a  circular  hole  (Figure  5). 

To  evaluate  the  performance  of  DSI  relative  to  the  MRG  method,  synthetic  experiments  were 
developed  by  applying  both  methods,  and  their  error  relative  to  the  known  analytic  solution  was 
calculated.  The  error  metric  utilized  was  defined  as 


el  =  le*  -s'  \) 
cm  I  °a  cm 

4  =  i4-4iJ’ 


(37) 


where  i  =  1 ...  N  is  indexing  of  the  N  evaluation  points,  ela  is  the  strain  as  calculated  from  the 
analytic  solution  at  point  i ,  and  s'm  and  s',  are  the  strain  on  point  i  from  the  MRG  and  DSI  methods, 
respectively.  The  mean  absolute  error  over  a  set  of  nodes  is  calculated  by 
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Figure  5.  Orthotropic  infinite  plate  with  circular  hole  and  load  at  infinity.  Direction  of  hatch  lines  coincides 
with  the  coordinate  system’s  Ox  axis  and  the  major  orthotropic  direction. 


em  -  n  iz 

i  =  1 

1  N 

i=  1 

Furthermore,  because  the  prevailing  source  of  errors  in  full-field  measurements  is  that  of  the  noise 
by  the  imaging  sensor,  a  simple  sensitivity  study  was  conducted  for  various  levels  of  noise.  To  inves¬ 
tigate  the  effect,  the  actual  strain  field  might  have  on  both  methods;  the  synthetic  tests  with  varying 
load  inclination  were  also  conducted. 


(38) 


6.1.  Orthotropic  plate  with  open  hole 

We  consider  an  infinitely  large  orthotropic  elastic  plate  that  contains  a  circular  hole  of  radius  r  that 
is  subjected  to  a  stress  state  at  infinity  equal  to  p.  The  applied  load  forms  an  angle  /3  with  the 
major  orthotropic  axis  as  shown  in  Figure  5.  This  angle  essentially  controls  the  degree  of  anisotropy 
introduced  in  the  strain-stress  fields  and  varying  it  amounts  to  rotating  the  major  orthotropic  axis 
relative  to  the  loading  direction  (i.e.,  dual  description  of  the  angle).  Strain  at  any  point  on  the  plane 
is  expressed  by  the  elastic  constitutive  equation,  [34] 

£xx  =  al\Oxx  +  0-12  Gyy 
£yy  —  & \2®xx  T  d22Qyy 
Yxy  —  @16®. xx  ~t~  a2(,Oyy 


+  a\6Txy  I 

+  Cl26Txy  i  ■  (39) 

+  0(,(,Txy  J 


It  has  been  shown  [34]  that  the  composition  of  Equation  (39)  for  the  general  case  of  multiply 
connected  plane  problems  can  be  reduced  to  a  set  of  algebraic  equations  in  terms  of  holomorphic 
functions  expressed  in  the  complex  plane  with  the  help  of  the  mathematical  theory  of  elasticity.  In 
fact,  the  general  solution  to  this  problem  for  an  elliptic  hole  is  known  [35]  and  is  implemented  here 
for  an  ellipse  of  minor  to  major  axes  aspect  ratio  of  unity,  to  reflect  the  problem  when  the  hole  is 
circular.  Additional  modifications  have  been  introduced  [36]  to  account  for  rotations  of  the  medium 
at  infinity. 
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The  stresses  are  calculated  by  equation  [35] 
Oxx  =  +  2Re 


230o(Ci)  ,  2dfo(Z2) 
Sl  - b52 


dzi 


dz  2 


(Jyy  =  a™  +  2Re 


TXy  —  TXy  2Re 


d<po{Zi)  dfo(Z2) 

dz\  dz.2 

dcpo (Zi)  ,  3iAo(Z2) 


Si- 


dzi 


+  S2~ 


dz.2 


(40) 


where  .v  i  and  S2  are  the  roots  of 

QllS^  —  2rti6 5^  +  (2fli2  “b  ^66)^~  —  2.Q26S  “b  ^22  =  0. 


(41) 


It  has  been  shown  [34]  that  Equation  (41)  has  no  real  roots  and  are  therefore  always  of  the  form 

.sT i  3  =  a\  ±  b\i,  52,4  =  a2  ±  bzi,  b\  >  0,  b2  >  0.  (42) 


The  stress  state  at  infinity  is  given  by 

axx  —  pco&2p,  a™  =  psin2/3,  r™  —  p  cos  /3  sin  P, 
whereas  the  holomorphic  function  is 


<Po(Zi)  = 

fo(z2)  = 


ir  [q_g°  +  is2o™  +  (s2  +  i)x™\ 
2(5i  -  s2) 

ir  [a™  +  is\o™  +  (5i  +  i) r™] 
2(51  ~s2) 


£t(~i) 

?2  (Z2) 


where 


(43) 


(44) 


Si(zi)  = 


&{Z2)  — 


Zi±  y/zi-r2(l+s2) 
r{ 1  +  /  5 1 ) 

C2  ±  ^z\~r2{\  +sj) 
r(l  +  is2) 


. 


(45) 


The  sign  ambiguity  in  Equation  (45)  is  removed  by  requiring  |£i|  $  1  and  |£2|  $  1.  The  com¬ 
plex  variables  are  defined  in  terms  of  the  real  position  coordinates  according  to:  £1  =  x  +  ,v  1  y 
and  z.2  —  x  +  52y  ■  Finally,  the  displacement  field  components  at  Ox  and  Oy  directions  are  given 
by  [36] 

u  =  2Re[ppp0(zi)  +  P2fo(z2)]  +  xe™  +  ye™  \ 

> ,  (46) 

v  =  2Re[qi<p0 (£1)  +  42^0(22)]  +  ye™  +  xe™  j 

with 


Pi  —  a\\S\2  +  cl\2  — ciif,Si 
P2  =  a11^22  +  Ol2  —  <216^2 


Cl\2S \2  +  0.22  ~  026^1 


0\2$22  +  O22  ~  026^2 


>  . 


(47) 
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Figure  6.  Strain  and  absolute  error  contours  for  nodal  noise  level  of  1  x  10— 2  pixels  and  loading  angle  of 
/I  =  45°.  The  domain  is  assumed  to  span  1600  pixels  wide  by  1600  pixels  in  height,  (a)-(c)  exx,  £yy,  £Xy 
of  the  analytic  solution;  (d)-(f)  £xx,  £yy,  £xy  of  the  synthetic  experiment  using  the  meshless  random  grid 
method;  (g)-(i)  £xx,  £yy,  £xy  of  the  synthetic  experiment  using  the  DSI  method;  (j)-(l)  respective  absolute 
error  of  the  meshless  random  grid  method;  and  (m)-(o)  respective  absolute  error  of  the  direct  strain  imaging 
method.  MRG,  meshless  random  grid;  DSI,  direct  strain  imaging. 
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6.2.  Numerical  experiments 

To  compute  the  fields  for  the  synthetic  numerical  experiments,  the  radius  of  the  hole  was  selected  to 
be  a  —  10 ~2m,  the  stress  at  infinity  a0 0  —lx  10 5 Pa,  the  angle  of  the  major  orthotropic  axis  was 
P  —  45°,  whereas  the  extent  of  the  domain  of  interest  (in  both  the  x  and  y  coordinates)  was  8a.  The 
material  properties  were  selected  to  be 


an 

am 

0 16 

2.245 

-2.732 

0.000 

am 

022 

026 

= 

-2.732 

9.436 

0.000 

a\6 

O22 

O  66 

0.000 

0.000 

4.315 

(48) 


A  total  of  2380  randomly  distributed  nodes  were  generated  using  the  method  described  in  [13], 
whereas  the  displacement  field  of  Equation  (46)  was  used  to  generate  the  displaced  nodes,  with 
varying  levels  of  noise.  Because  for  the  current  digital  imaging  technology  and  for  most  practical 
applications,  the  accuracy  in  coordinate  measurement  is  approximately  at  the  range  of  5  x  1 0~4  to 
5  x  10-2  pixels  [13];  the  noise  level  was  chosen  to  vary  from  10-4  to  10-1  pixels,  assuming  the 
entire  domain  is  mapped  on  a  digital  image  that  is  1600  pixels  wide  by  1600  pixels  high.  For  the 
load  angle  sensitivity  trials,  the  angle  /3  was  chosen  to  vary  from  0°  to  90°,  the  noise  level  was  kept 
at  1  x  10-2  pixels,  and  the  synthetic  runs  were  executed  70  times,  to  also  obtain  an  estimation  of  the 
variation  of  the  mean  absolute  error. 

The  basis  function  for  MRG  was  chosen  to  be  of  second  order  (a  total  of  six  terms),  whereas  for 
the  DSI,  it  was  chosen  to  be  of  first  order  (a  total  of  three  terms).  This  choice  was  made  so  that 
the  approximation  on  strains  for  both  methods  is  of  the  same  polynomial  degree.  The  domain  of 
support  radius  was  chosen  at  3.6  x  1  0~3m,  whereas  the  weight  function  was  the  one  presented  in 
Equation  (21). 

In  Figure  6,  the  strain  contours  as  well  as  the  absolute  error  of  the  MRG  and  DSI  methods  are 
plotted.  For  those  plots,  the  level  of  noise  imposed  on  the  node  coordinates  was  1  x  10-2  pixels.  It 
is  evident  that  DSI  performs  considerably  better  within  the  entire  domain  (Figures  6(m),  6(n),  and 
6(o))  compared  with  MRG  for  all  strain  components.  This  fact  is  more  evident  in  Figure  7,  where 
the  mean  absolute  error  of  exx  over  the  entire  domain  is  plotted  for  MRG  and  DSI  for  various  noise 
levels.  The  ratio  of  those  errors  (bottom  graph  of  Figure  7)  reveals  that  DSI  is  about  at  least  1.2 
times  more  accurate  than  MRG.  In  the  same  graph,  it  can  be  seen  that  DSI  can  be  up  to  3.3  times 
more  accurate. 

Similarly,  in  Figure  8,  the  mean  absolute  error  over  the  boundary  of  the  hole  is  plotted.  In  this 
case  as  well  DSI  outperforms  MRG  by  a  factor  ranging  from  about  1.2  to  about  3.6  times.  It  is  also 
worthwhile  mentioning  that  DSI  performs  better  with  increasing  noise  levels. 


Noise  level  [pixels] 

Figure  7.  Mean  absolute  error  of  exx  over  the  interest  domain  and  relative  performance  of  direct  strain 
imaging  over  meshless  random  grid  {eyy,  exy  show  similar  behavior). 
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Noise  level  [pixels] 


Figure  8.  Mean  absolute  error  of  exx  over  the  hole  boundary  and  relative  performance  of  direct  strain 
imaging  over  meshless  random  grid  ( syy ,  exy  show  similar  behavior). 


Figure  9.  Mean  absolute  error  of  exx  over  the  entire  domain  versus  load  angle  (syy  and  exy  show  similar 
behavior).  For  each  load  angle,  70  trials  were  performed,  and  their  standard  deviation  is  shown  as  the  error 
bars  in  the  graph.  It  is  worthwhile  noting  that  the  standard  deviation  of  the  direct  strain  imaging  error  is  very 
small  and  the  error  bars  are  almost  indistinguishable  from  the  trend  line. 


Both  methods  seem  to  reach  a  plateau  for  noise  levels  below  2  x  10-3,  which  is  anticipated  to 
depend  on  the  spatial  density  of  the  nodes.  It  should  be  noted  that  for  most  practical  cases  the  dis¬ 
placement  noise  levels  are  between  5  x  10-3  to  10  x  10-3  that  give  the  DSI  an  average  performance 
markup  of  about  3  times  for  both  error  metrics. 

Finally,  in  Figure  9,  the  mean  absolute  error  of  each  of  the  methods  is  plotted  with  respect  to  the 
loading  angle  as  this  induces  the  strain  variability  due  to  the  specific  nature  of  the  anisotropic  elastic 
problem  at  hand.  As  expected,  the  average  of  the  mean  error  is  not  affected  by  the  load  angle,  but 
for  the  MRG  method,  its  deviation  (i.e.,  the  error  of  the  error)  is  considerably  larger.  This  observa¬ 
tion  further  supports  the  utilization  of  the  DSI  method  in  applications  where  the  knowledge  of  the 
uncertainty  of  the  full-field  measurements  is  of  importance. 

7.  CONCLUSIONS 

In  this  work,  the  numerical  foundation  of  the  DSI  method  was  presented  and  its  accuracy  relative 
to  the  MRG  method  was  assessed.  DSI  outperformed  MRG  both  on  the  regions  near  the  domain 
boundaries,  where  full  field  methods  traditionally  suffer,  but  also  on  the  rest  of  the  domain  of  inter¬ 
est.  For  the  chosen  medium  configuration,  DSI  was  found  to  be  at  least  1.2  times  more  accurate, 
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while  there  were  tests  where  it  achieved  an  accuracy  improvement  of  up  to  3.6  times.  It  was  estab¬ 
lished  that  at  the  noise  levels  of  the  most  popular  digital  imaging  devices,  DSI  outperforms  MRG  by 
a  factor  of  3  in  terms  of  the  mean  absolute  error.  DSI  also  proved  to  be  very  robust  and  numerically 
stable  in  terms  of  mean  error  deviation,  as  was  shown  by  repeated  numerical  experiments  under  the 
same  image  sensor  noise  levels. 

Furthermore,  DSI  does  not  make  any  explicit  or  implicit  assumptions  on  the  continuity  of  the 
underlying  field,  and  therefore,  it  can  be  used  on  media  that  either  cannot  be  considered  continua  or 
have  developed  cracks  or  otherwise  failed. 

In  a  more  general  setting,  the  analysis  presented  herein  introduces  a  method  for  approximating 
tensor  quantities  with  respect  to  its  transformation  law,  from  directional  component  measurements, 
rather  than  a  scalar  quantity  that  is  usually  the  case  for  typical  regression  analysis  approaches. 

Future  plans  include  the  development  of  a  more  detailed  analysis  for  DSI  with  respect  to  the 
parameters  that  affect  its  performance  (noise  level,  node  density,  choice  of  polynomial  basis,  and 
choice  of  DOS),  as  well  as  comparison  of  its  performance  and  validation  with  actual  experimental 
data  obtained  by  strain  gauge  measurements. 
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